clc
clear
close all

thisFile = mfilename('fullpath');
thisDir = fileparts(thisFile);
cd(fileparts(fileparts(thisDir)));

load('data/processed/OilData.mat')

%% Table 1  Summary table: [min, mean, max, std. dev. ] 
% arrange the countries alphabetically
Q = OilData.Production(:,sort(OilData.Production.Properties.VariableNames));
% Price is in 2019:Q4 dollars per barrel
% Output is in 1000 barrel per day so we need to convert it to barrels per month by muliplying it by 30 * 1000 

Summary = zeros(4,20); %min, mean, max, std dev
format short
Output = zeros(size(Q));
for id=1:20
    q = Q{:,id}; % q is production in 1000 barrels per day
    q = q.*30.*1000./1000000; % in millions barrel per month
    Summary(1,id)=min(q);
    Summary(2,id)=mean(q);
    Summary(3,id)=std(q);
    Summary(4,id)=max(q);
    Output(:,id)=q;
    clear q
end
names = Q.Properties.VariableNames;
Summary     = array2table(Summary, 'VariableNames', names);
Summary.Price= [min(OilData.Price);mean(OilData.Price); std(OilData.Price); ...
    max(OilData.Price)];
Summary_min_mean_std_max = Summary;

% Save the table
writetable(Summary_min_mean_std_max, 'output/tables_and_figures/table1_summ_stats.csv');

%% Figure 1: time series of production   
for i=1:20 
    subplot(5,4,i)
    plot(OilData.Date, Output(:,i), 'LineWidth', 3)
    grid on
    set(gca, 'LineWidth', 1.5)
    ylim([0,max(Output(:))]) %creates same y-axis range 
    yticks([100,200,300,400,500])
    title(names{i})
end
saveas(gcf,'output/tables_and_figures/production_timeseries','epsc')
saveas(gcf,'output/tables_and_figures/production_timeseries','fig')
close all 
%% Figure 2: time series of prices
subplot(1,2,1)
plot(OilData.Date, OilData.Price, 'LineWidth', 3);
grid on 
set(gca, 'LineWidth', 1.5)
xlabel('Year')
ylabel('Price of Crude Oil/barrel')

subplot(1,2,2)
hold on 
histogram(OilData.Price,30, 'Normalization', 'pdf');
[t, xx]=ksdensity(OilData.Price,OilData.Price, 'Function', 'pdf', 'support', 'positive'); 
plot(xx, t,'LineWidth', 3) 
hold off 
grid on 
set(gca, 'LineWidth', 1.5)
xlabel('Price of Crude Oil/barrel')
ylabel('Probability')

saveas(gcf,'output/tables_and_figures/price_timeseries','epsc')
saveas(gcf,'output/tables_and_figures/price_timeseries','fig')
close all 
%% Figure 4 of Detrended Output by group
% detrending the data (x must be row vector)
T = size(Output,1); I = size(Output,2);
tdata = repmat((1:T)',1,I);
tempfun = @(x,tdata) -exp(repmat(x(1:I),T,1)-x(2*I+1)*tdata)...
    + exp(repmat(x(I+1:2*I),T,1));
tempfun1 = @(x) sum(sum((Output - tempfun(x,tdata)).^2));
options = optimoptions('ga','MaxGenerations',250*(2*I+1));
tau_est = ga(tempfun1,2*I+1,[],[],[],[],zeros(1,2*I+1),[9^999*ones(1,2*I) 1],[],[],options);

detrend_Q = Output + exp(repmat(tau_est(1:I),T,1)-tau_est(2*I+1)*tdata) ;

% GROUPING: 
% create group (for figure) in alphabetical order to match names
group = cell(6,1); 
group{1} = [7:9, 14, 17]; % iran, iraq, kuwait, qatar, uae
group{2} = [3,4,13,18]; % canada, china, norway, uk
group{3} = [5,11,20]; % Ecuador, mexico, venezuela
group{4} =[1,2,6,10,12];  %algeria, angola, egypt, libya, nigeria 
group{5} = [15,16]; % russia, saudi
group{6} = 19; %usa

for groupid=1:6
    subplot(2,3,groupid)
    hold on
    for country=1:numel(group{groupid})
        plot(OilData.Date,...
            detrend_Q(:,group{groupid}(country)),'LineWidth', 2)
    end
    title(['Group',num2str(groupid)])
    legend(names{group{groupid}},'Location','best')
    hold off
    grid on
    set(gca, 'LineWidth', 1.5)
    xlabel('Date')
    ylabel('(detrended) Output')
end

saveas(gcf,'output/tables_and_figures/detrendedQ','epsc')
saveas(gcf,'output/tables_and_figures/detrendedQ','fig')
close all 
